Entangled spin-orbital phases in the bilayer Kugel-Khomskii model 
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We derive the Kugel-Khomskii spin-orbital model for a bilayer and investigate its phase diagram 
depending on Hund's exchange Jh and the e g orbital splitting E z . In the (classical) mean-field 
approach with on-site spin (Sf) and orbital (rf ) order parameters and factorized spin-and-orbital 
degrees of freedom, we demonstrate a competition between the phases with either G-type or A-type 
antiferromagnetic (AF) or ferromagnetic long-range order. Next we develop a Bethe-Peierls- Weiss 
method with a Lanczos exact diagonalization of a cube coupled to its neighbors in ab planes by the 
mean-field terms — this approach captures quantum fluctuations on the bonds which decide about 
the nature of disordered phases in the highly frustrated regime near the orbital degeneracy. We 
show that the long-range spin order is unstable in a large part of the phase diagram which contains 
then six phases, including also the valence-bond phase with interlayer spin singlets stabilized by 
holes in 3z 2 — r 2 orbitals (VBz phase), a disordered plaquette valence-bond (PVB) phase and a 
crossover phase between the VBz and the A-type AF phase. When on-site spin-orbital coupling is 
also included by the {S*t* ) order parameter, we discover in addition two entangled spin-disordered 
phases which compete with A-type AF phase and another crossover phase in between the G-AF 
phase with occupied x 2 — y 2 orbitals and the PVB phase. Thus, the present bilayer model provides 
an interesting example of spin-orbital entanglement which generates novel disordered phases. We 
analyze the order parameters in all phases and identify situations where spin-orbital entanglement is 
crucial and mean-field factorization of the spin and orbital degrees of freedom leads to qualitatively 
incorrect results. We point out that spin-orbital entanglement may play a role in a bilayer fluoride 
K3CU2F7 which is an experimental realization of the VBz phase. 
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I. INTRODUCTION 

Recent interest and progress in the theory of spin- 
orbital superexchange models was triggered by the obser- 
vation that orbital degeneracy drastically increases quan- 
tum fluctuations which may suppress long-range order in 
the regime of strong competition between different types 
of ordered states near the quantum critical point P The 
simplest model of this type is the Kugel-Khomskii (d 9 ) 
model introduced long agc0 for KCuF 3 , a strongly cor- 
related system with a single hole within degenerate e g 
orbitals at each Cu 2+ ion. Kugel and Khomskii showed 
that many-body effects could then give rise to orbital or- 
der stabilized by a purely electronic superexchange mech- 
anism. A similar situation occurs in a number of com- 
pounds with active orbital degrees of freedom, where 
strong on-site Coulomb interactions localize electrons (or 
holes) and give rise to spin-orbital superexchangeP^The 
orbital superexchange may stabilize the orbital order by 
itself, but in e g systems it is usually helped by the orbital 
interactions whic h follow from the Jahn- Teller distortions 
of the latticePEHU For instance, in LaMn0 3 these contri- 
butions are of equal importance and both of them are 
necessary to explain the observed high temperature of 
the structural transition.^ Also in KCuF 3 the lattice dis- 
tortions play an important role and explain its strongly 



anisotropic magnetic and optical propertiesPHni 

An important feature of spin-orbital superexchange, 
which arises in transition metal oxides with active or- 
bital degrees of freedompnS is generic frustration of the 
orbital part of the superexchange. It follows from the di- 
rectional nature of orbital interactions, 1 which is in con- 
trast to the SU(2) symmetry of spin interactions. There- 
fore, the orbital part of the spin-orbital superexchange is 
intrinsically frustrated also on lattices without geometri- 
cal frustration, such as the three-dimensional (3D) per- 
ovskite lattice of KCuF 3 or LaMn0 3 . Generic features 
of this direction-dependent orbital interactions are best 
captured within the two-dimensional (2D) quantum com- 
pass modelpS which exhibits a quantum phase transition 
from one to the other one-dimensional (ID) columnar or- 
der through a poi nt with isotropic and strongly frustrated 
interactions! 12 * 14 * In spite of the intrinsic frustration and 
high degeneracy of the ground state, the long-range or- 
der of ID type exists in the 2D quantum compass model, 
as shown by a rigorous proofP^ Numerical simulations 
demonstrate that this model is in the universality class 
of the 2D Ising modePS and the order persists in a range 
of finite temperatureP21 In contrast, the superexchange 
interactions for the 2D e g orbital mode l contain orbital 
quantum fluctuations on the bondsp^ but nevertheless 
the long-range order survives also in this caseP^ 
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An intriguing situation arises when spin and orbital 
part of the superexchange are strongly coupled and com- 
pete with each other, as found in realistic spin -orbital 
models for several transition metal oxidesP^ For in- 
stance, a qualitatively new spin-orbital liquid phase may 
arise when the superexchange interactions are geomet- 
rically frustrated on the triangular lattice)^ or spin or- 
der cannot stabilize in LiNi02, another compound with 
triangular lattice of magnetic, in spite of presence of 
strong orbital interactions which suggest pronounced or- 
bital order A more standard situation is found in the 
transition metal oxides which crystallize in the perovskite 
lattice, where in general spin order coexists with or- 
bital order J^tSl and both satisfy the classical Goodenough- 
Kanamori rules PS A well known example is the archetyp- 
ical compound with degenerate e g orbitals, KCUF3, in 
which the orbital order is stabilized jointly b y th e su- 
perexchange and Jahn- Teller lattice distortions!^^ As a 
result, the magnetic interactions are strongly anisotropic 
and give rise to quasi-lD Heisenberg antiferromagnetic 
(AF) chain dominated by quantum fluctuations and char- 
acterized by spinon excitations, 21 with a dimensional 
crossover occurring when temperature is lowered below 
the Neel temperature T/yP^ 

While the coexisting A-type AF (A-AF) order and 
the orbital order is well established in KCuF 3 below 
T/vp^l and this phase is reproduced by the spin-orbital 
d 9 superexchange model)**' the model poses an interest- 
ing question by itself: Which types of coexisting spin 
and orbital order (or disorder) are possible when its mi- 
croscopic parameters are varied? So far, it was only es- 
tablished that the long-rang e AF order is destroyed by 
strong quantum fluctuations ) 24 ' 25 ' and it has been shown 
that instead certain spin disordered phases with valence- 
bond (VB) correlations stabilized by local orbital corre- 
lations are favoredP^ However, the phase diagram of the 
Kugel-Khomskii d 9 model is unknown — it was not stud- 
ied systematically beyond the mean-field (MF) approxi- 
mation and certain simple variational wave functions and 
it remains an outstanding problem in the theoryP 

The purpose of this paper is to analyze a simpler situ- 
ation of the spin-orbital Kugel-Khomskii model for a bi- 
layer, called below bilayer spin-orbital d 9 model, consist- 
ing of two ab layers connected by interlayer bonds along 
the c axis. This choice is motivated by an expected com- 
petition of the long-range AF order with VB-like states. 
One of them, a VB phase with spin singlets on the in- 
terlayer bonds (VBz phase), is stabilized by large crys- 
tal field E z which favors occupied 3z 2 — r 2 orbitals (by 
holes) . We shall investigate the range of stability of this 
and other phases, including the A-AF phase similar to 
the one found in KCuF 3 . 

To establish reliable results concerning short-range or- 
der in the crossover regime between phases with long- 
range AF or FM order, we developed a cluster MF ap- 
proach which goes beyond the single site MF in the spin- 
orbital systerrPS and is based on an exact diagonalization 
of an eight-site cubic cluster coupled to its neighbors by 



MF terms. This unit is sufficient for investigating both 
AF phases with four sublattices and VB states, with spin 
singlets either along the c axis or within the ab planes. 
This theoretical method is motivated by possible spin- 
orbital entanglement 27 which is particularly pronounced 
in the ID SU(4) [or SU(2)<g>SU(2)] spin-orbital models,^ 
and occurs also in the models for perovskites with AF 
spin correlations on the bonds where it violates the 
Goodenough-Kanamori rulesPHl In the perovskite vana- 
dates such entangled states play an important role in 
their optical properties,^ in the phase diagram^ and in 
the dimerization of FM in teract ions along the c axis in 
the C-AF phase of YV0 3 .£imi Below we shall investi- 
gate whether entangled states could play a role in the 
present Kugel-Khomskii model for a bilayer with nearly 
degenerate e g orbitals. Thereby we establish exotic type 
of spin-orbital order stabilized by joint quantum spin- 
orbital fluctuations, and investigate signatures of entan- 
gled states in this phase. 



The paper is organized as follows. In Sec.|TT]we present 
the Kugel-Khomskii d 9 spin-orbital model for a bilayer 
which consists of two 2D square lattices in ab planes cou- 
pled by vertical bonds along the c axis. First in Sec. |II A| 
we introduce the d 9 spin-orbital model for a bilayer de- 
rived here following Ref. 24, Its classical phase diagram 
obtained in a single-site MF approximation is presented 
in Sec. |IIB| Next we argue that quantum fluctuations 
and intrinsic frustration of the superexchange near the 
orbital degeneracy motivate the solution of this model in 
a better MF approximation based on an embedded cu- 
bic cluster, which we introduce in Sec. |II C| It leads to 
MF equations which were solved self-consistently in an 
iterative way, as described in Sec. IIP In Sec. |III| we 



present two phase diagrams obtained from the MF anal- 
ysis using Bethe-Peierls- Weiss cluster method: (i) the 
phase diagram which follows from factorization of spin 
and orbital degrees of freedom in Sec. Ill A and (ii) the 



one obtained when also on-site joint on-site spin-orbital 
order parameter is introduced, see Sec. |IIIB| The lat- 
ter approach gives nine different phases, and we describe 
characteristic features of their order parameters in Sec. 
|IV| We introduce bond correlation functions in Sec. |V A| 
and concentrate their analysis on the regime of almost 
degenerate e g orbitals, focusing on the proximity of the 
plaquette VB (PVB) and entangled spin-orbital (ESO) 
phases in Sees. |VB| and |V C| Finally, we quantify the 
spin-orbital entanglement using on-site and bond correla- 



tions, see Sec. VI which modifies significantly the phase 
diagram of the model with respect to the one obtained 
when spin and orbital operators are disentangled. Gen- 



eral discussion and summary are presented in Sec. VII 
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II. SPIN-ORBITAL MODEL AND METHODS 

A. Kugel-Khomskii model for a bilayer 

For realistic parameters the late transition metal ox- 
ides or fluorides are strongly correlated and electrons 
localize in the 3d orbitals j- 33 * 34 ! leading to Cu 2+ ions 
with spin S = 1/2 in d 9 configuration, as e.g. in 
KCuF 3 or La 2 Cu04. The virtual charge excitations lead 
then to superexchange which involves also orbital de- 
grees of freedom in systems with partly filled degenerate 
orbitals. In an alogy to the models introduced for bi- 
layer manganitej 35 ! 36 ! La2_ x Sr :! ;Mn207, we consider here 
a model for K3CU2F7 bilayer compound, with two active 
and nearly degenerate e g orbitals, 



(* 2 



V 



'■)/V2, \z)^(3z 2 -r 2 )/V6, (2.1) 



while t2 g orbitals do not contribute and are filled with 
electrons. They do not couple to e ff 's by hopping through 
fluorine and hence can be neglected. We investigate in 
what follows an electronic model and neglect coupling 
to the lattice distortions arising due to Jahn- Teller ef- 
fect. The bilayer K3CU2F7 system is known since twenty 
years but its magnetic properties were reported only 
recently.^ We shall address the orbital order and mag- 
netic correlations realized in this system below. 

The Hamiltonian for d 9 systems contains: holes' ki- 
netic energy H t with hopping amplitude t, electron- 
electron interactions i?; n t, with on-site Hubbard U and 
Hund's exchange coupling J H , as well as crystal-field 
splitting term H z playing a role of external orbital field 
E z acting on e g orbitals: 



H e „ = Ht + H\ 



int. 



H z 



(2.2) 



Because of the shape of the two e g orbitals Eq. 1 2.1 1, the 



effective hopping elements are direction dependent and 
different depending on the direction of the bond (ij). 
The only non-vanishing (dda) hopping element in the c 
direction connects two \z) orbitals, 6 while the elements 
in the ab planes satisfy Slater-Koster relations. 

Taking the effective (dda) hopping element t for two z 
orbitals on a bond along the c axis as a unit, H t is given 
by 



Ht = 



^ ] ^ d ixtr d jxa d ixa d jza) 



(ij) \\ab 



±V3(dL<W + 4xA za ) + H.c.} 
+ * E (dlAza + H-c), (2.3) 

fe> \\c 

where d\ xa and d\ za are creation operators for a hole in 
x and z orbital with spin a =t,4, and the in-plane x—z 
hopping depends on the phase of \x) orbital involved in 
the hopping process along the bond (ij ) and is included 
in the alternating sign of the terms cx v3 between a and 



b cubic axes. The on-site electron-electron interactions 
are described byP^ 



-Hint = U riiafTHai + (U - 3Jh) n ixa n lza 

ia ia 

+ (U- 2J H ) E n lxa n lzS - J H ^d\ x<J d ix3 d\ zS d lza 

ia ia 

+ JHj2( d U d U d ^i d ^ +H - c ^- ( 2 - 4 ) 

i 

Here rii aa stands for the hole density operator in orbital 
a = x, z with spin a, and a — — a . This Hamiltonian 
describes the multiplet structure of d s or d 2 ions and is 
rotationally invariant in the orbital space. We assumed 
the wave function to be real which gives the same am- 
plitude Jh for Hund's exchange interaction and for pair 
hopping term between \x) and \z) orbitals. The last term 
of the H Sg Hamiltonian lifts the degeneracy of the two 
e g orbitals 



H z ^ E z (jljxu Tliza) 3 



(2.5) 



and favors hole occupancy of x (z) orbitals when E z > 
(E z < 0). It can be associated with a uniaxial pressure 
along the c axis, induced by the bilayer geometry or by 
external pressure. 

The typical energies for the Coulomb U and Hund's 
exchange Jh elements can be deduced from the atomic 
spectra or found using density functional theory with 
constrained electron densities. Earlier studies performed 
within the local density approximation (LDA) gave 
rather large values of the interaction parameters!^! JJ = 
8.96 eV and Jh = 1.19 eV. More recent studies used the 
LDA with on-site Coulomb interaction treated within the 
LDA+J7 scheme and gave somewhat reduced values®' 
U = 7.5 eV and Jh = 0.9 eV. However, both param- 
eter sets give rather similar values of Hund's exchange 
parameter, 



n 



Jh 
U ' 



(2.6) 



being close to 0.13 or 0.12, i.e., within the expected range 
0.1 < 77 < 0.2 for strongly correlated late transition metal 
oxides. Note that the physically acceptable range which 
follows from Eq. (2.4) is much broader, i.e., < 77 < 1/3. 

The value of effective intersite (dda) hopping element 
t is more difficult to estimate. It follows from the usual 
effective process via the oxygen orbitals described by a 
t p d hopping, and the energy difference between the 3<i 
and 2p orbitals involved in the hopping process, so-called 
charge-transfer energyP A representative value of t ~ 
0.65 eV may be derived from the realistic parameters^ 
of Cu02 planes in La2CuC-4. Taking in addition U = 7.5 
eV, one finds the superexchange constant between hole 
5 = 1/2 spins within \x) orbitals in a single Cu02 plane, 
J x = (9/A)t 2 /U ~ 0.127 eV, which reproduces well the 
experimental value, as discussed in Ref. [231 
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Thanks to t <C U we can safely assume that the ground 
state is insulating at the filling of one hole localized at 
each Cu 2+ ion. In the atomic limit (t = and E z = 0) 
we have large 4 JV -fold degeneracy as the hole can occupy 
either x or z orbital and have up or down spin. This high 
degeneracy is lifted due to effective superexchange inter- 
actions between spins and orbitals at nearest neighbor Cu 
ions i and j which act along the bond (ij). They orig- 
inate from the virtual transitions to the excited states, 
i.e., d 9 d 9 ^dj°dj, and are generated by the hopping 
term Eq. (2.3). Hence, the effective spin-orbital model 



configuration, respectively, 
4 



S 7 ; • S. 



n 



m _ 



St -Si 0.10) 



for spins S = 1/2 at both sites i and j, and rf (with 7 = 
a, b, c standing for a direction in the real space) represent 
e g orbital degrees of freedom and can be expressed in 
terms of Pauli matrices {erf, a\ , erf} in the following way: 



a, b 



= -(-af±V3crf), 



(2.11) 



can be derived from the atomic limit Hamiltonian con- 
taining interaction Eq. (2.4) and the crystal-field term 
Eq. (2.5), treating the kinetic term Eq. (2.3) as a per- 



The matrices {erf} act in the orbital space (and have 
nothing to do with the physical spin Si present in this 
problem). Note that t? operators are not independent 



turbation. Taking into account the full mult ipletst rue- becauge they satigfy th * ^ constraint) £ 



r-'t = 



ture of the excited states for the d? configuration,^ 
gets the corrections of the order of Jh to the Hamilto- 
nian which results for the degenerate excited states (at 
Jh = 0). Calculating the energies of the excited d 8 states 
we neglected their dependence on the crystal-field split- 
ting E z . This assumption is well justified as the deviation 
from the equidistant spectrum at E z = become signif- 
icant only for \E z \/Jjj > 1 and in case of La2Cu04 one 
finds \E z \/Jh ~ 0.27. For systems close to orbital de- 
generacy, which we are interested in, this ratio is even 
smaller. 

The derivation which follows Ref. [pleads to the spin- 
orbital model, with the Heisenberg Hamiltonian for the 
spins coupled to the orbital problem, as follows: 



0. 



U = 



<*7>ll7 

+ (ra + r 4 ) 11™ 







*)] 


> 


(2.7) 



Here 7 = a, 6, c labels the direction of a bond (ij) in the 
bilayer system. The energy scale is given by the superex- 
change constant, 



J 



4*2 
IT 



(2.8) 



and the orbital operators at site i are given by 73 — 
{rf , t\ , rf} . The terms proportional to the coefficients 
{77, r 2 , t^} refer to the charge excitations to the upper 
Hubbard bancP^ which occur in the d 9 d 9 ^ d 9 d]j 9 pro- 
cesses and depend on Hund's exchange parameter 77 Eq. 
(2.6) via the coefficients!^ 



1 



ri 



I-377' 



1 — 7/ 



1 + 77 



(2.9) 



The model Eq. (2.7) depends thus on two parameters: 
(i) Hund's exchange coupling 77 Eq. (2.6), and (ii) the 
crystal-field splitting E z /J. 

The operators ILfj and stand for projections of spin 
states on the bond (ij) on a singlet (n^ ) and triplet (n*j ) 



In Fig. [T] we present typical orbitals configurations 
with ferro-orbital (FO) order and alternating orbi tal 
(AO) order considered in the e g orbital models In 
the next sections we shall analyze their possible coexis- 
tence with spin order in the bilayer d 9 spin-orbital model 
Eq. (2.7). As we can see, the maximal (minimal) value 
of the orbital operators rf is related with orbital taking 
shape of a clover (cigar) with symmetry axis pointing 
along the direction 7. 



B. Single-site mean-field approximation 

The bilayer spin-orbital d 9 model Eq. ( |2.7[ ) poses a 
difficult many-body problem which cannot be solved ex- 
actly. The only simple limits are either \E Z \ — ¥ 00 or 
77 — > (l/3)~ which we discuss below. In the first case 
the dominant term is the crystal field oc E z and, de- 
pending on its sign, we get uniform orbital configuration 
rf = ±1/2 and rf' b = =Fl/4. After inserting these classi- 
cal expectation values into the Hamiltonian Eq. \2.7\ we 
are left with the spin part which has purely Heisenberg 
form. 

We will show below that in the bilayer geometry of the 
lattice the single-site MF approximation predicts long- 
range ordered G-AF phases at 77 = known from the 
3D spin-orbital d 9 modelP see Fig. |2'd). For negative 
E z —> ~oo and FO order of z orbitals shown in Fig. [ljc), 
we get an AF coupling in the c direction and a weaker 
AF coupling in the ab planar directions (in the regime 
of small 77). For positive E z — > 00 one finds instead the 
FO order of x orbitals shown in Fig. [TJd), and two ab 
planes decouple, so we are left with the AF Heisenberg 
model on two independent 2D square lattices. In this 
case and the spins exhibit either G-AF, see Fig. |2jd) or 
G-AF order (not shown). Ferromagnetism is obtained in 
the present model for any E z if 77 is sufficiently large, i.e., 
when the superexchange is dominated by terms propor- 
tional to 7*1 which favor formation of spin triplets on the 
bonds accompanied by AO order depicted in Figs. [l|a) 
andQb). 

In what follows we will show the simplest, single-site 
MF approximation of the Hamiltonian Eq. (2.7) and 
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FIG. 1: (Color online) Schematic view of four representative 
orbital configurations on a cubic cluster: (a) AO order with 
(t"^ ) = 1/2 changing from site to site and (rf) = —1/4, ob- 
tained for E z < 0, (b) AO order with (t"' 6 '} = —1/2 changing 
from site to site and (rf) = —1/4, obtained for E z > 0, (c) 
FO order with occupied z orbitals and (rf) = —1/2 (cigar- 
shaped orbitals), and (d) FO order with occupied x orbitals 
and (t, c ) = 1/2 (clover-shaped orbitals). 



FIG. 2: (Color online) Schematic view of four different spin 
ordered phases on a cube realized in the d 9 spin-orbital model 
Eq. \2.7\ : (a) A-AF configuration, (b) plaquette valence 
bond (PVB) singlet configuration, (c) VBz phase with sin- 
glets along the c axis, and (d) G-AF configuration. Arrows 
stand for up or down spins, oval (violet) frames indicate sin- 
glets. Spin disordered phases with singlets on certain bonds 
(b) and (c) are stabilized by particular orbital order, see Sec. 
III. 



the resulting phase diagram. The Hamiltonian, originally 
expressed in terms of bond operators, can be then written 
in a "single-site" form given below: 



2,7 



(2.12) 



with sum running over all sites and cubic axes 7 = a, b, c. 
Here we adopted a shorthand notation with i+7 meaning 
the nearest neighbor of site i in the direction 7. 
The quantities 

x 7 = l ^mt^n-n % 7 = a ' 6 }, (2.13) 



and 



( r 2 - 

|(ra 



r 4 )H7 
h r 4 )W s 



if 7 = 0,6 
if 7 = c 



(2.14) 



are parameters obtained by averaging over spin opera- 
tors. The coefficients 1/2 in the x 7 and £ 7 terms along 
the c axis follow from the bilayer geometry of the lattice. 
We assumed that the spin order, determining x 7 and £ 7 , 
depends only on the direction 7 and not on site i. This is 
sufficient to investigate the phases with either AF or FM 
long-range order. More precisely, these are spin-singlet 



and spin-triplet projectors II 7 



= lT (t) 



(2.10) 



L s{t) — iJ -i,i+7 Ec L S 

that are independent of i. As far as only a single site 
is concerned the spins cannot fluctuate at zero tempera- 
ture and the projectors can be replaced by their average 
values: 



(Si ■ Si4 



n? 



(S,: • S 



1+7/ 



(2.15) 



The values of the projectors depend on the assumed spin 
order. Here we consider four different spin configurations: 
(i) G-AF - antiferromagnet in all three directions shown 
in Fig. [2](d) , (ii) C-AF - antiferromagnet in the ab planes 
with FM correlations in the c direction (not shown), (hi) 
A-AF - AF phase with FM order in the ab planes and AF 
correlations in the c direction depicted in Fig. [2ja), (iv) 
FM phase (not shown). The numerical values of the spin 
projection operators in these phases are listed in Table[T] 
After fixing spins, the MF approximation involves the 
well-know decoupling for the orbital operators: 



r 7 r 7 ~ (r 7 )r 7 +7 + r 7 (r 7 +7 ) - (r 7 )(r 7 +7 ). (2.16) 



'i 'i+7 



The last step is to define sublattices for the orbitals. The 
most reasonable choice would be to assume AO order 
meaning that neighboring orbitals are always rotated by 
7r/2 in the ab plane with respect to each other. To imple- 
ment this structure into the MF Hamiltonian we define 
new direction 7 as follows: 7 = 6, a for 7 = a, b and 7 = c 



() 



for 7 = c. Using 7 we can now easily define staggered 
order parameters: 



0.30 



tl 



i±7/ 



(2.17) 



The final single-site MF Hamiltonian can be written in 
the same form for any site so further on we will not use 
site index i anymore. The desired formula is: 

^Sf = E Q7r7 + = a ° z + + *<o> 



with 



and 



e 7 = ^ 7 +f(x 7 -P)-SAc, 



(2.18) 
(2.19) 



(2.20) 

For convenience we set J = 1; note that the energy scale 
can easily be recovered by replacing E z by E z / J. As we 
can see the MF Hamiltonian is very simple and can be 
written in terms of two Pauli matrices {a x ,a z } with 



1 



1 



1 



V3 



(2-21) 

Solving the 2x2 eigen-problem we obtain self-consistency 
equations for the order parameters t a and t c : 



2A 



(2.22) 
(2.23) 



where A = \J o? + f3 2 and the ground state energy given 
by: 



E 



-A-/(O c ). 



(2.24) 



The solution of self-consistency equations is very ele- 
gant and entertaining so we are not going to present it 
here and recommend it to the reader as an exercise (the 
results can be next compared with those given in the 
Appendix). It turns out that all four phases considered 



TABLE I: Mean values of triplet II/ and singlet HJ projection 
operators Eqs. ( |2.15[ ) for a bond (ij) along the axis 7 in 
different phases with long-range magnetic order which occur 
in the MF phase diagram, see Fig. [3] 



phase 


n a(b) 


n? 


n a(6) 


n° 


G-AF 


1/2 


1/2 


1/2 


1/2 


G-AF 


1/2 


1 


1/2 





A-AF 


1 


1/2 





1/2 


FM 


1 


1 










FIG. 3: (Color online) Phase diagram of the bilayer d 9 spin- 
orbital model Eq. (2.181 obtained in the single-site MF ap- 
proximation with spin and orbital MFs. In this approach 
the G-AF and G-AF phases (for E z > —0.25 and moderate 
rj) have exactly the same energy. Shaded gray (green) area 
indicates phases with AO order, while the remaining states 
with long-range G-AF spin order are accompanied by FO 
order. The enlarged area around the multicritical point at 
E z — — 0.25J and 77 = is shown in the inset. 



here can appear as orbitally uniform, i.e., having FO or- 
der with orbitals being either perfect clovers or perfect 
cigars everywhere, or as phases with AO order between 
two sublattices. The phase diagram presented in Fig. [3] 
was obtained by purely energetic consideration and shows 
the boarder lines between phases with the lowest energies 
for given 77 and E z . This diagram is surprisingly complex 
taking into account the simplicity of the single-site ap- 
proach; it reveals seven different phases. For r\ — we 
have only two AF phases: (i) G-AFz for E z < —1/4 J 
and (ii) G-AF for E z > —1/4 J, with a different but uni- 
form orbital configuration (FO order) which involves ei- 
ther cigar-shaped z orbitals in the G-AFz phase, see Fig. 
T 
T 



c) , or clover-shaped x orbitals in the G-AF, see Fig. 

d) . Because of the planar orbital configuration in the 
latter G-AF phase one finds no interplane exchange cou- 
pling and thus this phase is degenerate with the G-AF 
one. 

For higher 77 the number of phases increases abruptly 
by three phases with AO configurations, as shown in 
the inset of Fig. [3} the A-AF, G-AF/AO and G- 
AF/AO phase. Surprisingly, the AO version of the G- 
AF phase is connected neither to z nor to x FO order 
in an antifcrromagnet, excluding the multicritical point 
at (E z /J,ij) = (—0.25,0), and disappears completely for 
X] « 0.118. The G-AF/AO phase stays on top of uni- 
form G(G)-AF phase, lifting the degeneracy of the above 
phases at relatively large r\ and then gets replaced by the 
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FM phase which always coexists with AO order. One can 
therefore conclude that the G(C)-AF degeneracy is most 
easily lifted by turning on the orbital alternation. 

On the opposite side of the diagram the G-AFz phase 
is completely surrounded by A-AF phases: for rj > 
(2/v3 — 1) the G-AFz phase turns into orbitally uni- 
form A-XFz independently of the value of E z (interor- 
bital triplet excitations dominate then on the bonds in 
the ab planes), and for smaller 77 into the A-AF phase 
with AO order. In the A-AF phase the AF correlations 
in the c direction survive despite the overall FM tendency 
when rj grows. This follows from the orbitals' elongation 
in the c direction present for E z < 0, which would cause 
intcrplane singlets formation if we were not working in 
single-site MF approximation, see Sec. III. In the present 
case it favors either the G-AFz or A-AF(z) configuration 
with uniform or alternating orbitals depending on the 
values of E z /J and rj. Finally, the FM phase is favorable 
for any E z if only r\ is sufficiently close to 1/3 which only 
confirms that the single-site MF approximation is sound 
and not totally wrong with this respect. 

The central part of the presented diagram is the most 
frustrated one judging by the number of competing 
phases with long-range spin order. This behavior is con- 
sistent with that found in the 3D spin-orbital d 9 model 
in the regime of E z ~ and finite 17 P Four of these 
phases could be expected by looking at the phase di- 
agram of the 3D model: two G-AF phases, the A-AF 
phase and the FM phased Note, however, that in the 
phases stable in the central part of the phase diagram, 
namely in the A-AF, A-AF/AO and FM phase, the oc- 
cupied orbitals alternate. While the FM phase is not 
surprising in this respect and obeys the Goodenough- 
Kanamori rule of having FM spin order accompanied by 
the AO order, in the A-AF one finds an example that 
both spin and orbital order could in principle alternate 
between the two ab planes. This finding suggests that 
in this central part of the phase diagram one may ex- 
pect either other VB-type phases or even states with 
more complex spin-orbital disorder. Such ordered or dis- 
ordered phases require a more sophisticated approach, 
either variational wave functions pEU or the embedded 
cluster approach which we explain below in Sec. 



into identical cubes which cover the lattice, the Hamilto- 



nian (2.7) can be written in a cluster MF form as follows, 



Hmf = 



mec 



(U'mt , -i/cxt\ 



(2.25) 



where the sum runs over the set of cubes C, with indi- 
vidual cube labeled by C m <E C. Here T-V^ contains all 
bonds from T-L m belonging to the cube C m and the crystal 
field terms oc E z , i.e., it depends only on the operators on 
the sites inside the cube, while "H™* contains all bonds 
outgoing from the cube m and connecting neighboring 
clusters, making them correlated. 

The basic idea of the cluster MF approach is to ap- 
proximate by containing only operators from 
the cube m. This can be accomplished in many different 
ways depending on which phase we wish to investigate. 
Our choice will be to take HT? of the following form: 



n/c 

ft. 



+ r? $ + (2.26) 



II C 



containing spin field Sf breaking SU(2) symmetry, or- 
bital field T 2 7 and spin-orbital field Sfrf. Coefficients 
{aj, tij, cj, cf-} are the Weiss fields and should be fixed 
self-consistently depending on E z and 77. Our motivation 
for such expression is simple: if orbital degrees of free- 
dom are fixed then the problem reduces to the Heisenberg 
model which has long-range ordered AF phase — that 
is why we take Sf field, the orbitals are present in the 
Hamiltonian so taking r 7 is the simplest way of treating 
them on equal footing to describe possible orbital order. 
Finally, we introduce also spin-orbital field SjtJ because 
we believe that in some phases spins and orbitals alone do 
not suffice to describe the symmetry breaking and these 
operators can act together. 

The standard way to go on is to write self-consistency 
equations for the Weiss fields. This can be done in a 
straightforward fashion: we take the operator products 
from and divide them into a part depending only on 
Oi operators from the cube m, and on Oj ones — from 
a neighboring cube n. Then we use the well known MF 
decoupling for such operator products, 



C. Cluster mean— field Hamiltonian 

Now we introduce a more sophisticated approach 
which goes beyond the single-site MF approximation of 
Sec. |IIB| In what follows we use a cluster MF approach 
with a cube depicted in Fig. |4j It contains eight sites cou- 
pled to its neighbors along the bonds in ab planes by the 
MF terms. This choice is motivated by the form of the 
Hamiltonian with different interactions along the bonds 
in three different directions — the cube is the smallest 
cluster which does not break the symmetry between the 
a and b axes and contains equal numbers of a, b and c 
bonds. After dividing the entire bilayer square lattice 



O.O, « (OJOj + OiiO^-iOJiOj) 

= O t {0 3 ) - ^OiXO;) + 0,(0,:) - \(Oi)(Oi) , 

(2.27) 

and write it in a symmetric way. Now the first two terms 
can be included into W™' and the last two into "H° xt . 
This procedure can be applied to all operator products 
in H^* an d full "H™* can be recovered in the form given 
by Eq. (2.261. Repeating this for all clusters leads to a 



Hamiltonian describing a set of commuting cubes inter- 
acting in a self-consistent way. After using Eq. (2.27) 



on the Hamiltonian Eq. (2.7) we obtain the formulas for 




7 

m,i ' 



FIG. 4: Schematic view of the cluster used in the Bethe- 
Peierls- Weiss MF approach of Sec. |II C| Vertices i = 1, • • • , 8 
and directions 7 = a,b,c are marked in the figure, and dashed 
lines stand for the outgoing MF interactions in the a and b 
direction. 



the Weiss fields 
1 



7 



I 



0*2 + r 4 )u/ + -(r 2 - n)s 
1 



b] = -(r 4 +ri)u7--(r 2 -ri)s7, 



(2.28) 
(2.29) 
(2.30) 



7+7 

raA 



1 , 

- ^(3ri+2r 2 +r 4 ), 
where the order parameters at site i are: 



t 1 = 

u m,i — 



(S?), 



S : 



(2.31) 



(2.32) 
(2.33) 

(2.34) 



7 u 1 



} are the mean values of operators 



Note that {s^, t 
at site i belonging to the cluster m, and {sj ,ij ,u]} are 
the mean values of the same operators at sites neighbor- 
ing with i in the direction 7. The geometry of a bilayer 
implies that each site i has one neighbor along the axis a 
and another one along the axis b, and these sites belong 
to different cubes. 

The next crucial step is to impose a condition that 
{s],t],u]} are related to the order parameters obtained 
on the internal sites of the considered cluster. The sim- 
plest solution is to assume that all clusters have identical 
orbital configuration; t] — spin configuration is in 
agreement with a type of global magnetic order we want 
to impose; sj = ±Si and spin orbital configuration is 



as if spin and orbitals were factorized, i.e., uj = ±u 
This solution has one disadvantage: if a or b direction is 
favored in the orbital configuration of the cube then this 
broken symmetry will propagate through whole lattice 
which is contradictory with the form of the Hamiltonian 
Eq. ( |2.7[ ). That is why it is better to assume that two 
neighboring cubes can differ in orbital (and spin-orbital) 
configuration by the interchange of a and b direction, i.e., 



±Si 



tl 



(2.35) 



with 7 being the complementary direction in the ab plane 
to 7, i.e., (7,7) = (a, £>),(&, a). This relation gives the 
same results as the previous one in case when the (a, b) 
symmetry in the cube is not broken, but keeps the whole 
system (a, b) symmetric in the other case. Here we again 
treat the spin-orbital field as factorized but surprisingly 
it turns out that this does not prevent spin-orbital entan- 
glement to occur, see below. We have also tried to impose 
relations between u] and i which have nothing to do 
with spin and orbital sectors alone but this only resulted 
in the lack of convergence of self-consistency iterations. 



D. Self— consistent iterative procedure 

The self-consistency equations cannot be solved ex- 
actly because the effective cluster Hilbert space is too 
large even if we use total S z conservation in the consid- 
ered cluster m (then the largest subspace dimension is 
d = 17920) and because of their non-linearity. The way 
out is to use Bethe-Peierls- Weiss method, i.e., set cer- 
tain initial values for the order parameters {«i,C,w^J 
and next employ Lanczos algorithm to diagonalize "Hmf 
Eq. (2.25). Below we present results obtained by self- 



consistent calculations of phases with broken symmetry 
or with spin disorder. In order to determine the ground 
state one recalculates mean values of spin, orbital and 
spin-orbital fields given by Eqs. ( |2.35 ) and determines 
new order parameters. This procedure is continued until 
convergence (of energy and order parameters) is reached. 
This process can be very slow due to the number of or- 
der parameters which is 24 (three per site) for the cube 
but we have overcome this problem by imposing certain 
symmetry breaking on the cluster. Wc implement it 
in the following way: after each iteration we calculate 
{s^ 8 -} only for one site i = 1 and the remaining 

coefficients are fixed assuming certain symmetries of the 
phase we are searching for. 

For simplicity let us enumerate the vertices i = 1, ■ • • , 8 
in the cubic cluster as shown in Fig. [4] To obtain G-AF 
phases we assume that: 



si if i € A 
-si if i € B 



(2.36) 



for a two-sublattice structure, where A — {1,4,5,8} 
and B = {2,3,6,7}. In FM case it is enough to put 
Si = si and in case of FM order within the planes and 
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AF between them (in the A-AF phase) we use instead: 
Si = (— l) 4_1 si. In the orbital sector we can impose a 
completely uniform configuration with t'J ni = t^ nl , which 
can however lead to non-uniform configuration of the 
whole system because neighboring clusters are rotated 
by 7r/2 with respect to each other, or we can produce a 
phase with AO order taking: 



t! 



C,i if i^B 



(2.37) 



with (7, 7) = (a, b), (b, a). Other choices would be to take 
the above equation either with A = {1, 2, 5, 6} and B = 
{3,4,7,8} or with A = {1,3,5,7} and B = {2,4,6,8}. 
More generally speaking, every choice of orbital sublat- 
tices is good as long as the total MF wave function does 
not violate the symmetry between directions a and b. 
The sublattices for spin-orbital field are constructed as if 
1 could be expressed as i = - -f m i ). 
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III. PHASE DIAGRAMS 

A. Disentangled spin and orbital operators 

The zero-temperature phase diagr am of the present bi- 
layer d 9 spin-orbital model Eq. (2.7) depends on parame- 



FIG. 5: (Color online) Phase diagram of the d 9 spin-orbital 



ters {E z , rf) , and was obtained by comparing ground state 
energies for different sublattices formed by {s,, t , ik^ A- 
MFs. In this way we determined the ground state with 
the lowest energy and its order parameters. We begin 
with the phase diagram of Fig. [5] obtained by assuming 
that spin-orbital operators may be factorized into spin 
and orbital parts, i.e., i = Sj(l/2 — J or: 



{SI 



(3.1) 



Next we report the phase diagram (in Sec. IIIBl, where 
we include i calculated following the definition in Eq. 
(2.34). Comparing these two schemes allows us to de- 



termine which phases cannot exist without spin-orbital 
entanglement. 

The low-77 part of the diagram in Fig. [5] is dominated 
by three phases: VBz for negative E z , PVB for E z close 
to zero and G-AF for positive E z . The VBz phase with 
ordered interlayer valence bonds for occupied z orbitals 
and spin singlets, see Fig. [2^c) , has replaced the G-AFz 



phase obtained before in Sec. II B Both phases exhibit 
uniform FO order, i.e., i is close to —1/2 for all i which 
means that orbitals take the shape of cigars aligned along 
the c bonds, see Fig. [ijc). One finds that quantum 
fluctuations which could be included within the present 
approach select the VBz phase and magnetization van- 
ishes due to the singlets formation. For higher values 
of E z ~ also a different phase is found: the plaquette 
VB (PVB) phase with singlets formed on the bonds in a 
or b direction of the cluster, see Fig. ||b). This phase 
breaks the a—b symmetry of the model locally but the 
global symmetry is preserved thanks to the ir/2 rotation 



superexchange model Eq. (2.25) obtained using the cluster 
method for an embedded cube with factorized spin-orbital 
operators. Valence-bond phases with spin disorder are stable 
in the shaded area. 



of neighboring clusters (see Eq. (2.35)). The orbitals are 
again uniform within the cluster with tj^ i ovt b mi close to 
— 1/2, meaning that they take shape of cigars pointing in 
the direction of the singlets. For high positive values of 
E z the ground state is the G-AF phase with long-range 
AF order and FO order of x occupied orbitals, i.e., t c mi 
close to 1/2, see Fig. [ljd). This means that orbitals are 
indeed of the x type and take shape of four-leaf clovers in 
the ab plane with lobes pointing along a and b directions 
which makes the two planes very weakly coupled. 

The FO order in the VBz and G-AF phases agrees with 
the limiting configurations for E z — > ±00 described ear- 
lier. The first of them is a quantum phase with local sin- 
glets, in contrast to the G-AFz one found already in Sec. 



II B 



and in the 3D spin-orbital d 9 modelf^ If we consider 
now the VBz phase and increase 77, we pass through the 
VBm phase (where "m" stands for mixed orbital config- 
uration) and reach the ^4-AF phase with non-zero global 
magnetization such that spins order ferromagnetically in 
the ab planes and antiferromagnetically between them 
(along c axis), see Fig. [2]Ja). We believe that this regime 
of the phase diagram is of relevance for the spin and or- 
bital correlations in K3CU2F7 and discuss it also in Sec. 
VII[ The orbital order is of the AO type with t c m i close to 
zero, positive or negative depending on E z , see Figs, [l^a) 
andQb). The VBm phase occurs when the orbitals in 
the VBz phase start to deviate from the uniform configu- 
ration and ends when the global magnetization appears, 
accompanied by the change of the orbital order. The first 
transition is of second order, being the only second order 
phase transition in this diagram of Fig. [5j 

The presence of both A-AF phases on top of the VBz 
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can be understood qualitatively as follows: in the VBz 
phase AF spin coupling is strong only within the singlets, 
so if rj is increased the weak in-plane spin correlations can 
easily switch to FM ones, while AF correlations will still 
survive between the planes. The last phase of the dia- 
gram is the FM phase with AO order, similar to the AO 
order in the A-AY phase. Due to the absence of ther- 
mal and quantum fluctuations the magnetization in this 
phase is constant and maximal. The FM phase appears 
for any E z , if only rj is sufficiently close to 1/3, which 
agrees qualitatively with the previous discussion of the 
exact limiting configurations and with the phase diagram 
found before in the single-site MF approach, see Fig. [3] 
Comparing Fig. [5] to the MF phase diagram of Fig. 
[3] we can immediately recognize the main difference: the 
existence of the VBz and PVB phases. These phases con- 
tain spin singlets on the bonds and do not follow from the 
single-site MF approach. Another difference is the lack of 
sharp transitions between AO and FO order within one 
phase; these transitions are smoothened by spin fluctua- 
tions absent in the single-site MF and perfect FO config- 
urations are now available only for extremely high values 



B. Phase diagram with spin-orbital field 
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FIG. 6: (Color onli ne) T he phase diagram of the cluster MF 
Hamiltonian Eq. (2.251 of the d 9 spin-orbital model for a 



bilayer, with independent spin, orbital and spin-orbital mean 
fields. Valence-bond phases with spin disorder are stable in 
the light shaded (yellow) area, and phases with spin-orbital 
entanglement are indicated by dark gray (orange) shading. 



When the spin-orbital MF is not factorized but cal- 
culated according to its definition given in Eq. (2.341, 
one finds the phase diagram displayed in Fig. [6] We 
would like to emphasize that this non-factorizability can- 
not be included within the single-site MF approach be- 
cause there all spin fluctuations are absent. Of course, 
one can imagine that we take the Sf(Sf + ^} decoupling in 
the pure-spin sector and Sf(Sf +J ) decoupling in the spin- 
orbital sector of the Hamiltonian Eq. ( |2.7| ) leading to the 
fluctuating spins but this would break both the magne- 
tization conservation and homogeneity of the spin-spin 
interactions included into the Kugel-Khomskii model. 

In addition to the phases obtained in the phase dia- 
gram of Fig. [5j we get here also the following phases: 
ESO, EPVB and PVB-AF (the VBm phase is still stable 
between the VBz and A-AY ones but has much smaller 
area) . The first two above phases are formed in the highly 
frustrated region of the phase diagram where both E z 
and rj are moderate. ESO stands for entangled spin- 
orbital phase and is characterized by relatively high val- 
ues of spin-orbital order parameters, especially for high 
r\ values when other order parameters are close to zero. 
This phase contains singlets along the bonds parallel to 
the c axis, its magnetization vanishes and the orbital 
configuration is uniform. One can say that this is the 
VBz phase with weakened orbital order transformed into 
uniform spin-orbital order for the same spin and orbital 
sublattices. EPVB stands for entangled PVB phase and 
resembles it, but has in addition finite non-uniform spin- 
orbital fields, and weak global AF order. A different type 
of phase with spin-orbital entanglement is the PVB-AF 



phase connecting PVB and G-AF in a smooth (as it will 
be shown below) way but only if r\ is large enough. In 
contrast to the direct PVB-to-G-AF transition, passing 
through the PVB-AF involves second order phase transi- 
tions and the same happens in case of the EPVB connect- 
ing the ESO and PVB phases. Similarly to the previous 
diagram, the transition from the VBz to VBm phase is 
of the second order while the other transitions produce 
discontinuities in order parameters (see Sec. IV) and cor- 
relation functions (see Sec. [Vb. 

Finally, we should also point out that the G-AY/C- 
AF degeneracy found in Fig. [3] is lifted in the cluster 
approach and the C-AF phase does not appear in any 
of the two phase diagrams presented in Figs. [5] and [6] 
Another interesting feature of the phase diagrams are 
points of high degeneracy where different phases have the 
same ground-state energies. In case of the single-site MF 
diagram this quantum critical point is found at (E z — 
— 1/4 J, rj = 0), where six phases meet. The use of cluster 
MF method which includes singlet phases lifts this point 
upwards along the border line between VBz and PVB 
to (E z ,rj) « (—0.3 J, 0.11) in case of Fig. [6j This means 
that singlet formation acts against interaction frustration 
caused by Hund's exchange coupling and moves the most 
frustrated region of phase diagram to high-?7 regime. This 
shows once again that the simple single-site approach is 
not sufficient to describe correctly the properties of the 
bilayer d 9 spin-orbital model. 
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FIG. 7: (Color online) Order parameters {s, t a ' b , v a ' b } for 77 = 
0.05 and -0.4 < E z /J < 0.3 in the VBz, PVB and G-AF 
phase, from left to right. 



IV. THE ORDER PARAMETERS 

The ground state is characterized by order parameters 
obtained directly during the self-consistency iterations in 
each phase: spin, orbital and spin-orbital order parame- 
ters, {si, t c ^' 1 , w^j il- We focus here on the phases shown 
in the phase diagram of Fig. [6j For the physical reasons 
it is however better justified to define joint spin-orbital 
order parameter in a slightly different way, introducing a 



new variable v m i as follows: 



>.T — 



(4.1) 



which differs from the old order parameter by a subtrac- 
tion of the spin field, i.e., i = \si — v'l ni . Now one can 
study the behavior of order parameters along different 
cuts of the phase diagram of Fig. [6] and determine types 
of phase transitions. Below we present a few represen- 
tative results. For this purpose we first choose 77 = 0.05 
and start within the VBz phase, where by increasing E z 
one gets first into the PVB and next to G-AF phase, see 
Fig. [7} For 77 = 0.15 there are even more phases and 
one passes through the A-AF, ESO, EPVB, PVB, and 
AF-PVB phases, before reaching finally the G-AF phase, 
see Figs. [8] and [9| We also investigated the dependence 
of order parameters on Hund's exchange coupling — we 
fixed E z = — 0.72J, started in the VBz phase and in- 
creased 77 to get to the VBm and A-AF phases — these 
results are shown in Figs. [T0| 

In what follows we use shorthand notation for the order 
parameters, 



{S,t >,«'} = {«!,* ' i,« 'i> 



(4.2) 



In Fig. [7] we displayed the order parameters for increas- 
ing E z in phases VBz, PVB and G-AF (separated by 
dotted lines in the plot). The sublattice magnetization 



s is non-zero only in the G-AF phase because the re- 
maining phases are of the VB crystal type, with spin 
singlets oriented either along the c direction or in the ab 
planes. In the G-AF phase the spin order grows stronger 
for increasing E z when the orbital fluctuations weaken 
and spin fluctuations present in the G-AF phase reduce 
s from the classical value of 1/2. 

Consider now decreasing values of E z in Fig. [7j Both 
orbital order parameters remain equal and close to — 1/4 
in the G-AF phase until the (first order) transition point 
to the PVB phase, where orbital configuration changes 
abruptly and becomes anisotropic. In this case the a—b 
symmetry was broken in such a way that that spin sin- 
glets point in the PVB phase in b direction and so the 
directional orbitals (cigars) do. This explains the robust 
orbital order with t b being close to —1/2 in most of the 
PVB phase. The global symmetry is not broken as the 
VB singlets form here a checkerboard pattern in the ab 
plane, with AO order of directional orbitals along the a 
and b axis in the neighboring plaquettes. The transition 
to the VBz phase is discontinuous (first order) in the or- 
bital sector too: t a grows constantly while decreasing E z 
down to 0.4J, drops slightly close to the transition point 
and jumps to 1/4 in the VBz phase, t b grows quickly 
to t b pa 0.125 while approaching the transition and then 
jumps to the value of t a . Qualitatively this means that 
close to the above transition the orbital cigars pointing 
along the b axis change gradually into a shape very sim- 
ilar to clover orbitals lying in the be plane and then sud- 
denly the lobes along the b direction disappear and we 
are left with the pure VBz phase. 

The spin-orbital order parameter behaves in a much 
less intriguing way; it remains zero in the VBz and PVB 
phases, jumps to finite value at the PVB-to-G-AF tran- 
sition and remains almost constant and close to —0.1 in 
the G-AF phase. The vanishing value of v a,b in the sin- 
glet phases is simple to understand: the orbitals are here 
fixed and spins form singlets and fluctuate independently 
between the values ±1/2. This means that t a,b and s are 
not "synchronized" in any way and only this could lead 
to v a ' b ^ 0. This condition is satisfied in the G-AF phase; 
orbitals are fixed and the spin configuration is here de- 
termined by s order parameter. 

Figure [8] shows that the transition between the PVB 
and G-AF phases can have a completely different char- 
acter than described above. The difference comes from 
a higher value of rj which is now equal to 0.15, enhanc- 
ing the FM channel of superexchange and leading to the 
intermediate PVB-AF phase and to a smooth transition 
from the PVB to G-AF phase. In the PVB-AF phase 
staggered magnetization s grows continuously from zero 
(in the PVB) to a finite value in the G-AF phase and re- 
mains saturated there. This means that planar singlets 
in the PVB phase decay gradually and spins get par- 
tially "synchronized" with orbitals, moving toward uni- 
form configuration which gives finite spin-orbital order 
parameters v 1 ^ 0. The anisotropy (v a ^ v b ) follows 
from the anisotropy of orbitals inherited from the PVB 
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FIG. 8: (Color online) Order parameters {s, t a ' b , v a - b } for 77 = 
0.15 and 0.3 < E-JJ < 0.5 in the PVB, PVB-AF and G-AF 
phase, from left to right. 
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FIG. 10: (Color online) Order parameters {s,t a 



>} for 



E z = -0.72J and 0.15 < E z /J < 0.25 in the VBz, VBm and 
A-AF phases, from left to right. 




FIG. 9: (Color online) Order parameters {s, t a ' b , v a ' b } for 77 = 
0.15 and -0.5 < E z /J < -0.1 in the A-AF, ESO, EPVB and 
PVB phases, from left to right. 



order. This mechanism of the PVB-to-G-AF transition 
is absent for low values of 77 — we anticipate that the en- 
hanced FM component of interactions reduces spin fluc- 
tuations which makes the correlations between spins and 
orbit als possible. 

In Fig. [9]we focus on the complementary regime of the 
phase diagram, 77 = 0.15 and negative E z . In this regime 
we describe three different consecutive phase transitions 
between the phases: A-AF, ESO, EPVB and PVB. The 
first phase transition can be regarded as a little bit artifi- 
cial because this is a meeting point of two completely dif- 
ferent types of spin and orbital order, with different sym- 
metries and sublattices. For this reason the transition 
has to be discontinuous and the spin order parameter has 
different physical meaning on both sides of the transition 
line, i.e., s in the magnetic moment in the A-AF phase 
while it is a weak AF order parameter in the ESO phase. 



We anticipate that a smooth crossover occurs in place 
of such a transition in the thermodynamic limit, never- 
theless by comparing the energies we concluded that this 
transition follows from the cluster MF approach. Note 
also that the ESO phase has predominantly z orbitals 
accompanied by fluctuations, i.e., t c ~ —0.4 and t a = t b , 
and may be seen as an extension of the VBz phase. 

On the contrary, the second quantum EPVB phase 
which occurs in the phase diagram of Fig. [6] may be 
seen as a precursor of the PVB phase and is character- 
ized again by finite joint spin-orbital fluctuations, with 
v a =/= for a = a, b. What is especially peculiar in 
the EPVB phase is the non-zero staggered magnetiza- 
tion s which grows smoothly from the zero values at the 
phase borders meaning that we have a wedge of antiferro- 
magnetism between two VB configurations. The EPVB 
phase seems to be similar to PVB-AF in a sense that 
spin-orbital fields are non-zero and non-uniform but the 
qualitative behavior of the order parameters is different, 
e.g. in the EPVB phase spin-orbital fields have always 
opposite signs, while in the PVB-AF phase the sings are 
the same. 

Looking at the orbital order parameters t a,b in the A- 
AF phase (Fig. [9]), one observes similar anisotropy as in 
the PVB one but this time a-b symmetry is not broken 
within the cluster because in the A-AF phase every or- 
bital is rotated by 7r/2 with respect to its neighbors in the 
ab plane. Another difference is that the orbitals take the 
shape of clovers, not cigars, with symmetry axes pointing 
along the a or b axis which is described by t b being close 
to 1/2. In the A-AF phase we have also long-range mag- 
netic order and finite spin-orbital fields, indicating joint 
behavior of spin and orbital MF variables. 



Next Fig. 10 shows the behavior of order parameters 
for E z — —0.72 J and 77 changing in an interval allowing 
us to study the transitions from the VBz to VBm phase, 
and between the VBm and A-AF phase. In this case 
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FIG. 11: (Color online) Order parameters {s,t a ' b ' c ,v a ' b } for 
77 = 0.28 and -4.0 < E z /J < 6.0 in the A-AF, FM and G-AF 
phases, from left to right. 



all the phases can be described by the same spin and 
orbital sublattices because VBz is uniform in the orbital 
sector and has no long-range magnetic order so it can be 
described both in terms of the PVB and A-AF type of 
ordering. Global magnetization appears only in the A- 
AF phase jumping from the zero value in the VBm and 
growing with increasing 77. Transition from the VBz to 
VBm phase is continuous in both spin and orbital sectors. 



10 



The orbital order parameters t a ' h bifurcate in Fig. 
at r\ ~ 0.169 from the isotropic value t a = t b ~ 1/4 and 
the orbital anisotropy grows in the VBm phase to give 
AO order in the A-AF phase (Fig. 



10) 



and next shows 
a discontinuity at the second transition. The final AO 
order can be described by clover orbitals with symmetry 
axes alternating between a and b directions from site to 
site. Relatively big, negative value of t a means that the 
clovers' lobes are elongated in the a or b direction, per- 
pendicular to their axes. The elongation depends also 
on the value if E z : the E z — > —00 limit corresponds 
to pure clover-like orbitals, while for E z —> 00 one gets 
pure cigars. This tendency is especially visible in the 
FM phase which is not limited in horizontal direction of 
the phase diagram. Consequently, the VBm phase can 
be regarded as a crossover regime between orbitally uni- 
form VBz and alternating A-AF phases. This resembles 
to some extent the PVB-AF phase described earlier but 
we want to emphasize the main difference between these 
phases: the VBm phase does not need non-factorizable 
spin-orbital MF to appear while the PVB-AF one needs 
it (compare Figs. [5] and [6]). The question of spin-orbital 
non-factorizability will be addressed in more details be- 
low, see Sec. IVII 



Finally, we show the behavior of the order parame- 



ters {s,t 



a,b,c „,a.b 



} and the quantum fluctuation effects on 



on t a and t b (by the constraint t c — —t a — t b ), and was 
added here to visualize the orbital order along the c di- 
rection which is essential in the large \E Z \ regime showed 
in Fig. [TT] In FM phase the spin order is saturated be- 
cause of the lack of quantum and thermal fluctuation. For 
the same reason spin-orbital field factorizes and {v a ,v b } 
fields bring no extra information which would not be al- 
ready contained in {t a ,t b }. The overall behavior of t c is 
in agreement with the crystal field part of the Hamilto- 
nian Eq. (2.7) with t c —¥ ±1/2, giving uniform cigar or 



clover orbitals depending on the sign of E z . 

We emphasize that for increasing E z one finds two 
crossing points of t c with {t a ,t b } curves, one at t c = 
t b = -1/4 and the other one at t c = t a = 1/4. At 
these two points the orbitals take shapes of perfect clovers 
(E z < 0) or perfect cigars (E z > 0), with symmetry axes 
alternating in the ab plane from site to site. Only one 
of these points belongs to the FM phase meaning that 
the four "perfect" orbital configurations: AO order with 
clovers/cigars and FO order with clovers/cigars are sep- 
arated by phase transitions in the spin-orbital model Eq. 
(2.7). The transitions shown in Fig. 11 are discontinu- 
ous due to the change of global spin order in each phase. 
The spin order parameter s plays a role of staggered A- 
AF or AF magnetization in the extremal phases and is 
trivial (saturated) in the FM phase. On the other hand, 
all three phases displayed in Fig. 11 can be described 



them in the A-AF, FM and G-AF phases for 77 
Fig " 



11 The third orbital field t c 



0.28, see 
is linearly dependent 



by the same orbital sublattices assuming AO order. The 
large scale of E z in Fig. [TT|is in contrast to those in other 
figures — it indicates that orbital degrees of freedom are 
very rigid when spins are almost frozen and one needs 
rather high energies to change their configuration. 



V. NEAREST NEIGHBOR CORRELATIONS 



A. Spin, orbital and spin-orbital correlations 

Studying order parameters in different phases we get 
complete information about symmetry broken or disor- 
dered phases of the system, but this alone does not justify 
the use of the cluster MF method as order parameters 
can in principle be obtained using standard single-site 
MF approximation, see Sec. |IIB| The advantage of the 
cluster method becomes evident when we investigate cor- 
relation functions on the bonds belonging to the consid- 
ered cube. The most obvious ones are the spin-spin cor- 
relations (Si • Sj) or orbital-orbital correlations (t^tJ), 
but in addition one may also determine joint spin-orbital 
correlations, ((S^ • S^t^tJ). Although one could in prin- 
ciple invent several other bond correlation functions, the 
above ones have the most transparent physical meaning 
because they enter the Hamiltonian. For the same rea- 
son we will only consider orbital correlation functions for 
different bond direction 7 = a, b, c. This gives nine corre- 
lation functions, three in each direction, for each vertex 
of the cube. For symmetry reasons it is enough to con- 
sider only one chosen vertex, e.g. vertex number 1 in Fig. 
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|H For convenience we will use the following notation: 



CI = (Si -Si 

r<1 — 
u st = 



'1 

((Si 



Si)riV N 



(5.1) 
(5.2) 
(5.3) 



where the bond (li)||7 and i € {2,3,7} which gives all 
nonequivalent nearest neighbor correlations along 7 = 
a, 6, c (see Fig. |3|. 

In the next paragraphs we will present the numerical 
results for bond correlations {C] , C7, C] t } along different 
cuts of the phase diagram of Fig. [6j For all three-panel 
plots each panel describes correlation of one type: upper 
panel — spin correlations, middle one — orbital correla- 
tions, and bottom one — spin-orbital correlations. For 
each panel different characters (colors) of line indicate 
different direction 7: solid (red) line stands for 7 = c, 
dashed (green) line for 7 = and dashed-dotted (blue) 
line for 7 = 6. In case of two-panel plots there are only 
two directions considered, c and a because for symmetry 
reasons correlations along the b and a axes are identical. 
Therefore, the left panel concerns all three types of cor- 
relators for 7 = c and the right one for 7 = a in the way 
that solid (red) lines are spin-spin correlation functions, 
dashed (green) ones are orbital-orbital correlations and 
dashed-dotted (blue) represent spin-orbital correlators. 
In order to investigate the nature of spin-orbital entan- 
glement we focus the discussion on two quantum phases 
which occur at finite values of Hund's exchange rj near 
the orbital degeneracy: (i) the PVB phase, and (ii) the 
ESO phase. 



B. Plaquette valence-bond phase 

We begin with bond correlation functions for rj = 0.05 
and -0.4 < E z /J < 0.3 in the VBz, PVB and G-AF 
phase. The G^r function stays close to —3/4 in the VBz 
phase while the other spin correlations are almost zero 
as one can expect in the interlayer singlet phase, see Fig. 



12 a) . After the first transition at E z ~ — 0.26J the 
situation changes — now the singlets are in b direction 
and C b s gets close to —3/4 when E z increases. After the 
second transition at E z ~ 0.14 J all the spin correlations 
take finite negative values with C c s relatively weakest, 
keeping the symmetry between a and b direction. This 
is in agreement with the spin order in the G-AF phase 
discussed in Sec. \FV\ 

The orbital correlation functions in the VBz and G-AF 
phases behave as if the orbitals were frozen in uniform 
configuration with t c = ±1/2 and t a,b = ±1/4 whereas 
in the PVB phase their behavior is more nontrivial; the 
dominant C b is quite distant from its maximal value 1/4 
and the difference between G t a and G t c is visible, espe- 



cially close to the G-AF phase, see Fig. 12 Id). This result 
is due to quantum fluctuations: perfect VBz and G-AF 
configurations are the exact eigenstates of the Hamilto- 
nian, at least in the limit of large \E Z \, while perfect PVB 
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FIG. 12: (Color online) Nearest neighbor correlations for rj = 
0.05 and -0.4 < E z /J < 0.3 in the VBz, PVB and G-AF 
phases, (a) Spin correlations: solid (red) line — CJ, dashed 
(green) line — C" and dashed-dotted (blue) line — C*. (b) 
Orbital correlations: solid (red) line — Cf, dashed (green) 
line — Cf and dashed-dotted (blue) line — C\. (c) Spin- 
orbital correlations: solid (red) line — Cg t , dashed (green) 
line — Cg t an d dashed-dotted (blue) line — C st . 



state cannot be obtained exactly in any limit and gets 
easily destabilized by varying E z . It is peculiar that the 
spin configuration is almost nonsensitive to the orbital 
splitting E z and the singlets stay rigid in the regime of 
spin disordered phases, i.e., below the transition to the G- 
AF phase. 



The spin-orbital sector, shown in Fig. 12 'c), 
does not bring any new information; all the lines behave 
as if spin and orbital degrees of freedom were factorizable. 



Figure 13 presents the bond correlations for a grad- 
ual transition between the PVB and G-AF phases, with 
an intermediate PVB-AF phase for 77 = 0.15 and 0.3 < 
E z /J < 0.5. By decreasing E z , i.e., looking from right 
to left, we can see the in-plane spin correlation bifur- 
cating smoothly at the transition to the PVB-AF phase 
and evolving monotonically to the values characteristic 
of the PVB phase, see Fig. 13 'a). The interplane spin 



correlations G^ stay relatively weak everywhere which is 
obvious in both PVB and G-AF phase and hence not so 
surprising in the intermediate PVB-AF phase. 

In the orbital sector we can see here very similar be- 
havior to the one observed in Fig. 12 — again the order 
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FIG. 13: (Color online) Nearest neighbor correlations for 77 = 
0.15 and 0.3 < E z /J < 0.5 in the PVB, PVB-AF and G-AF 
phases, (a) Spin correlations: solid (red) line — CJ, dashed 
(green) line — C" and dashed-dotted (blue) line — C h a . (b) 
Orbital correlations: solid (red) line — C t c , dashed (green) 
line — Ci and dashed-dotted (blue) line — G\. (c) Spin- 
orbital correlations: solid (red) line — Cg t , dashed (green) 
line — and dashed-dotted (blue) line — C h st . 
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FIG. 14: (Color online) Nearest neighbor correlations for 77 = 
0.15 and -0.45 < E z /J < -0.1 in the ESO, EPVB and PVB 
phases, (a) Spin correlations: solid (red) line — CJ, dashed 
(green) line — C" and dashed-dotted (blue) line — Cg. (b) 
Orbital correlations: solid (red) line — Cf, dashed (green) 
line — Ct and dashed-dotted (blue) line — C\. (c) Spin- 
orbital correlations: solid (red) line — Cg t , dashed (green) 
line — Cst an d dashed-dotted (blue) line — C st . 



is far from the perfect PVB but C t c is close to the classical 
value of 1/16 obtained for the plane perpendicular to two 
directional orbitals along the b axis, while C t ° is almost 
exactly opposite and C\ stays below 1/4, see Fig. 13 1)). 



This shows some kind of universality at the transition 
from the PVB to G-AF phase which is independent of 
the intermediate phase. Again, the spin-orbital sectors, 
shown in Fig. 13 c) , does not indicate any qualitatively 
new behavior comparing to spins and orbitals alone but 
looking at the phase diagrams with (Fig. [5]) and with- 
out (Fig. [6]) spin-orbital factorization we recognize that 
on-site spin-orbital entanglement must be responsible for 
the onset of the PVB-AF phase. 



C. Phases with entangled spin-orbital order 

Consider now smaller (negative) values of E z , where 
unexpected and qualitatively new entangled phases occur 
in the phase diagram of Fig. [6j We display bond correla- 
tion functions in Fig. 14 in two neighboring highly frus- 



trated and entangled phases, the ESO and EPVB phase 
- the latter one turns into the PVB phase when E z is 
increased. The relevant parameter range for r\ = 0.15 
is -0.45 < E z /J < —0.1. On the first glance this plot 
shows that the transitions between the ESO and EPVB 
as well as between the EPVB and PVB phases are of the 
second order. In the spin sector one observes weaken- 
ing singlet order in the ESO phase with C c s getting far 
from —3/4 and in-plane correlations C^ h being practi- 
cally vanishing, see Fig. 14 '&). After the first transi- 
tion (at E z ~ —0.36 J) C b s grows rapidly toward negative 
values while C c s goes to zero much more gently and C° 
stays close to zero. This means that in the EPVB phase 
we have relatively strong AF order in the be plane inside 
the cluster, turning into the ac plane order on neigh- 
boring cubes. This gives finite magnetization s shown 
in Fig. [9j When approaching the second transition (at 
E z ~ -0.22J) C c s weakens and C b s gets closer to -3/4 
and this is continued within the PVB phase. 

In the orbital sector we can find other differences be- 
tween entangled and disentangled phases, see Fig. [Tlfcb). 
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FIG. 15: (Color online) Nearest neighbor correlations for 
E z = -0.43 J and 0.15 < r] < 0.25 in the ESO and A-AF 
phases, (a) Correlations along the c axis: solid (red) line — 
Cg, dashed (green) line — C£ and dashed-dotted (blue) line 

— Cst- (b) Correlations within ab planes: solid (red) line — 
Cg , dashed (green) line — Cf and dashed-dotted (blue) line 

— Cf f . Correlations in a and b direction are the same. 
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FIG. 16: (Color online) Nearest neighbor correlations for 
E z = -0.72J and 0.10 < rj < 0.25 in the VBz, VBm and A- 
AF phases, (a) Correlations along the c axis: solid (red) line 

— Cg, dashed (green) line — Ct and dashed-dotted (blue) 
line — Cg t . (b) Correlations within ab plane: solid (red) line 

— Cg , dashed (green) line — C t a and dashed-dotted (blue) 
line — Cg t - Correlations in a and b direction are the same. 



In the ESO phase the drops considerably when ap- 
proaching the first transition; this is in contrast with the 
VBz phase where C t c stays almost constant until the tran- 
sition occurs. However, one finds that the spin-orbital 
bond correlation C c st stays constant in the ESO phase, 
see Fig. 14 ^c). The behavior of in-plane correlation func- 



tions C t a ' becomes somewhat puzzling within the EPVB 
phase: after bifurcation at the transition point C\ drops 
to zero and slowly recovers to become dominant in the 
PVB phase, while C t a stays dominant in certain region 
of the EPVB phase even though the spin correlations 
in a direction vanish. Only C t c gradually drops to zero 
throughout all three phases. 

Note that in the spin-orbital sector we can see the joint 
order in both entangled phases in a more transparent way 
than in the orbital one, at least concerning the ESO and 
EPVB phases (we should keep in mind that —3/16 < 
C] t < 1/16 while -1/4 < C"l < 1/4 where the bottom 
limit for C] t is realized only in singlet phases). The C c st 
correlation is definitely dominant in the ESO phase and 
stays dominant in most of the EPVB phase in contrary 
to spin correlation. In addition, close to the second 
transition the C c st correlation is overcome by C\ t which 
grows here stronger because of singlets being formed on 
the bonds along the b axis. This tendency is further 
amplified within the PVB phase. Note that C° t stays 
practically zero in all the phases shown in Fig. [14] 

Now we turn to the dependence of bond correlations 
on increasing Hund's exchange rj. In Fig. [15] we display 
correlations for E z = — 0.43J and 0.15 < tj < 0.25 in the 
ESO and A-AF phases. Both phases can be described by 



a strong tendency toward AO order with two sublattices 
which does not violate the a-b symmetry inside the cube; 
for this reason we show only correlations along the c and a 
direction. The spin sector within the ESO phase is dom- 
inated by the decay of interplanar singlets accompanied 
by growth of in-plane correlations which triggers global 
A-AF order above the transition (at -q ~ 0.22). The or- 
bital correlations in the c direction drop almost to zero 
when rj grows and stay small in the A-AF phase. The in- 
plane orbital correlations C t a ' b decrease in the ESO phase 
too but remain finite after the transition. Summarizing, 
in the ESO phase close to the onset of the A-AF one we 
find a very weak orbital order accompanied by precursors 
of the A-AF order in spin sector. 

Consider now the spin-orbital correlations. In the ESO 
phase C c st takes relatively big, negative values and does 
not change much except for the transition point where 
it jumps to zero. In contrast, in the A-AF phase we 
no longer observe any spin-orbital ordering. Note that 
a peculiar signature of the ESO phase is rather robust 
spin-orbital order on the interlayer bonds along the c 
axis which turns out to be more rigid against quantum 
fluctuations than orbital order and remains finite even 
when orbital order vanishes. 

In the last figure, Fig. [l6j we display bond correlation 
functions in the VBz, VBm and A-AF phases for E z — 
-0.72J and 0.10 < rj < 0.25. As before, all the in- 
plane correlations are independent of 7. The plots prove 
that the transition from the VBz to VBm phase is of the 
second order while the transition from the VBm to A-AF 
phase produces no discontinuities in correlations either, 
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but the behavior of order parameters (see Fig. 10) is 
here slightly discontinuous. In the spin sector we observe 
first (at r\ < 0.17) that robust singlets along the c axis 
with ~ —0.7, see Fig. [l6^a), are gradually weakened 
under increasing rj and weak FM correlations occur in the 
VBz phase close to the first phase transition to the VBm 
order. We suggest that this regime of parameters could 
correspond to K3CU2F7, where the magnetic properties 
indicate interplanar singlets as formed in the VBz and 
VBm phases accompanied by weak FM correlations in 
the ab planes 

Note that the changes in spin correlations with increas- 
ing T] become fast only after leaving the VBm phase. In 
the orbital sector perfect VBz order dies out quickly al- 
ready in the VBm regime, both on the bonds along the c 
and a axes. After entering the A-AF phase, G t c vanishes 
exponentially while G t a crosses zero and gradually falls to 
negative values. This behavior is in agreement with that 
shown in Fig. [TO] saying that t c remains close to zero in 
the A-AF phase and the negative Gf confirms AO order 
in ab planes. Altogether, the spin-orbital sector does not 
exhibit here any considerable non-factorizable features. 



VI. SPIN-ORBITAL ENTANGLEMENT 

The essence of spin-orbital entanglement observed 
in the cluster MF approach is spin-orbital non- 
factorizability. This feature can have either on-site or 
bond character, the latter was introduced in Ref. [27] 
We emphasize that on-site entanglement which is charac- 
teristic for cases with finite spin-orbit coupling, 42 occurs 
also in the present superexchange model as shown below. 
We define the on-site entanglement as non-separability 
of the order parameters, i.e., spin and orbital operators 
are entangled when v 7 ^ st 7 while the entanglement 



as being of bond type wherP-0 C] t 7^ C]C^ , implying 
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that it can be detected by investigating the respective 
correlation functions. Therefore we analyze in this Sec- 
tion the numerical results for the quantities (covariances) 
motivated by the above discussion which are defined as 
follows: 



r 7 = v 1 — st 1 , 



(6.1) 
(6.2) 



In case of r 7 we consider only 7 = a, b as the on-site 
covariance satisfy the local constraint, 

r c = -r a - r b , (6.3) 

while for R 1 we shall present the data for 7 = a, b, c. 
In order to quantify the above non-factorizability and to 
recognize whether it is strong or weak in a given phase, 
it is necessary to establish first the minimal and max- 
imal values of R? and r 7 . Simple algebraic considera- 
tions give the following inequalities: the bond covariances 
|i? 7 | < 0.25 in singlet phases, |i? 7 | < 0.125 in phases with 
magnetic order, and the on-site covariances |r 7 | < 0.25 
everywhere. 

First of all, the numerical results show that both bond 



Eq. (6.2 ) and on-site Eq. (6.1 ) spin-orbital entanglement 



E z < 0.3 in the VBz, PVB and G-AF phases. 



is small in the regime of weak Hund's exchange coupling. 
This feature is illustrated in Fig. [17] for the VBz, PVB 
and G-AF phases at 77 = 0.05 and -0.4 < E z /J < 0.3. 
The r° = r b curves show no on-site spin-orbital entan- 
glement (r 1 — 0) in both VBz and PVB phases, while 
it is finite in the G-AF phase (r a = r b < 0) and grad- 
ually approaches zero with increasing E z . We empha- 
size that this on-site non-factorizability is minute, being 
one order of magnitude smaller than its maximal value, 
and does not play any important role for the stability of 
the G-AF ground state. This is confirmed by the fact 
that G-AF phase exists in the same region of parameters 
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in both phase diagrams: factorizable (Fig. [5]) and non- 
factorizable one (Fig. [6]), and occurs even in the single- 
site MF approximation (Fig. [3|. It is interesting to note 
that the in-plane bond entanglement R a ' b takes relatively 
high values in the G-AF phase. This is clearly an effect of 
quantum fluctuations; the perfect (classical) G-AF phase 
of Fig. [3] has uniform fixed x orbital configuration with 
t c = 1/2 which suppresses any non-factorizability. As the 
on-site entanglement, also the bond spin-orbital entan- 
glement vanishes gradually for high values of E z — > oo. 

At the border line between the VBz and PVB phases 
we noticed a considerable increase of R c and less pro- 
nounced growth of R b which seem to be induced by the 
transition as the R b,c drop quickly for higher values of 
E z . In the VBz phase we expect all the spin-orbital co- 
variances to be zero for the same reasons as in the G-AF 
phase and this also applies to the perfect PVB phase. In 
Fig. [TTJ however, the VBz and PVB phases are domi- 
nated by the critical behavior which distorts perfect or- 
der ings. 

Also in the regime of higher Hund's exchange inter- 
action 77 = 0.15 the spin-orbital covariances in the PVB, 
PVB-AF and G-AF phases are small in the range of their 
stability, see Fig. [l8]for 0.3 < E z /J < 0.5. In the PVB 
phase all the covariances take small values showing that 
the PVB type of order has no serious quantum fluctu- 
ations in this parameter range. The on-site covariances 
{r a ,r b } bifurcate from the zero value at the first transi- 
tion and this emergence of non-factorizability stabilizes 
here the intermediate PVB-AF phase (compare Figs. [5] 
and [6]) and persists in the G-AF phase where they over- 
lap again (r a — r b ). In the regime of PVB-AF phase 
we observe also almost linear decrease of the in-plane 
R a ' b staying close to each other and a smaller drop of 
R c . Although these quantities are all small, the order 
parameters (see Fig. |8| are small too, so we conclude 
that spin-orbital entanglement is qualitatively important 
here. The minimum of all R 1 is located at the second 
transition indicating that highly entangled states play a 
role also at the onset of the G-AF phase. 

Figure [19] shows spin-orbital entanglement in the most 
exotic part of the phase diagram with the ESO, EPVB 
and PVB phases for ?y = 0.15 and -0.45 < E z /J < -0.1. 
The on-site spin-orbital covariances {r a ,r b } take high, 
opposite values in both the ESO and EPVB phase, 
with maximum (minimum) at the transition line between 
them. Comparing to other phases r a b values are high- 
est in the ESO and EPVB phases, and comparing the 
two phase diagrams in Figs. [5] and [6j we recognize that 
spin-orbital entanglement is a constitutive feature of both 
ESO and EPVB states. We emphasize that the on-site 
spin-orbital entanglement is strong and complementary 
in the ESO phase on the bonds along the a and b direc- 
tion (r a = — r b ), while it vanishes between the ab planes 
(r c = 0). These results indicate spin-orbital fluctuations 
in the ab planes, with (S z a x ) ^ and no fluctuations 
along the c axis, where r c follows from (S z a z ) = 0. In 
contrast, in the EPVB phase there is also finite on-site 
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entanglement for the interlayer order parameters, r c 0. 

Looking at the bond parameters K 1 we see that the 
dominant one is R c falling gradually in the ESO down to 
the minimum at the ESO-EPVB transition. At the same 
point R b drops from zero value in the ESO phase and 
takes maximally negative value inside the EPVB regime. 
In contrast, R a remains close to zero in the entire regime 
of parameters and in the PVB all the covariances go to 
zero showing that the order within the PVB phase is 
practically disentangled. The dominant role of R c comes 
from the c-axial symmetry of the ESO phase and in- 
creased quantum fluctuations on the ESO-EPVB border 
while the non-zero value of R b in the EPVB phase follows 
from the magnetic and orbital order on the cube in the 
be plane mentioned in the previous section. 

When Hund's exchange is increased across the transi- 
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tion between the ESO and phases, one finds that 

bond and on-site spin-orbital covariances are radically 
different in both phases, see Fig. 20 The plot shows 



that the ESO phases is much stronger entangled than 
the ^4-AF one where all the covariances stay close to zero. 
Only the in-plane R a ' b parameters are small also in the 
ESO phase but the other covariances, including the bond 
covariance along the c axis R c , take considerable values. 

Finally, we focus on the range of large negative crystal 
field splitting E z — — 0.72J and display the spin-orbital 
covariances in the VBz, VBm and ^4-AF phases for in- 
creasing Hund's exchange 0.15 < r\ < 0.25, see Fig. [21] 
On the one hand, looking at the VBm region of the plot 
we can understand why this phase can exist when factor- 
ized spin-orbital MF is applied (again, compare Figs. [5] 
andpj); the on-site covariances {r a ,r b } vanish here and 
within the VBz phase. On the other hand, one finds cer- 
tain on-site entanglement in the ^4-AF phase, especially 
close to the transition line — this shows why the A-AF 
area is expanded in Fig. [6] as compared with the non- 
factorized phase diagram of Fig. [5] Concerning bond 
entanglement, it is significant (finite R c < 0) only along 
the interlayer c bonds in all these three phases, taking 
maximal values of |i? c | in the A-AF phase. One can un- 
derstand this as follows: in the VBz phase the orbital 
order is almost perfect and orbitals stay frozen — there- 
fore spin-orbital factorization is here almost exact as in- 
dicated by a low value of R c . This is not the case in the 
A-AF phase where orbitals fluctuate, especially close to 
the transition line to the VBm phase. The R a ' b bond pa- 
rameters are small due to the imposed FM order within 
the ab planes which decouples the spin from orbital fluc- 
tuations on the bonds along the a and b directions. 



VII. DISCUSSION AND SUMMARY 

The numerical results presented in the last three sec- 
tions, obtained using the sophisticated mean-field ap- 
proach with an embedded cubic cluster, provide a trans- 
parent and rather complete picture of possible ordered 
and disordered phases in the bilayer spin-orbital d 9 
model. This approach is well designed to determine the 
character of spin, orbital and spin-orbital order and cor- 
relations in all the discussed phases as it includes the 
most important quantum fluctuations on the bonds and 
captures the essential features of spin-orbital entangle- 
ment. By analyzing order parameters we also presented 
evidence which allowed us to identify essential features 
of different phases and to distinguish between first and 
second order phase transitions. This is especially impor- 
tant in cases when two phases are separated by an in- 
termediate configuration, such as the PVB-AF or VBm 
phase, where one finds a gradual evolution from singlet to 
AF correlations, supported (or not) by non-factorizable 
spin-orbital order parameter. We believe that the cluster 
mean-field approach presented here and including thie 
joint spin-orbital order parameter is more realistic be- 
cause there is no physical reason, apart from the form of 
the d 9 Hamiltonian, to treat spin and orbital operators 
as the only fundamental symmetry-breaking degrees of 
freedom in any phase. 

Interestingly, the results show that the bottom part 
of the phase diagram of the d 9 spin-orbital model does 
not exhibit any frustration or spin-orbital entanglement 
up to r\ ~ 0.07 and the type of spin and orbital order 
are chosen there predominantly by the crystal field term 
oc E z . Quantum fluctuations dominate at E z ~ and 
for E z < 0, where they stabilize cither in-plane or inter- 
planar spin singlets accompanied by directional orbitals 
which stabilize them. These two phases arc reminiscent 
of the resonating valence-bond (RVB) phase and PVB 
phase found in the 3D spin-orbital d 9 modelPHere, how- 
ever, the VBz phase extends down to large values of E z , 
where instead the long-range order in the G-AFz phase 
was found in the 3D model. This demonstrates that inter- 
layer quantum fluctuations are particularly strong in the 
present bilayer case. On the contrary, at E z > one finds 
the G-AF spin order which coexists with FO order of x 
orbitals. It is clear that both the VBz and G-AF phase 
are favored by the interplay of lattice geometry and by 
the shape of occupied orbitals for E z — > ±oo. In this low- 
•q regime of the diagram the area occupied by the PVB 
phase is narrow and especially orbital order is affected 
by the quantum critical fluctuations. The planar singlets 
are formed shortly after leaving the VBz phase and re- 
main stable afterwards. Spin-orbital non-factorizability 
seems to be marginal in the entire VBz phase but plays 
certain role when switching to the planar singlet phase, 
especially visible for the interplane bond covariance R c . 

On the contrary, in the PVB phase away from the crit- 
ical regime spin-orbital non-factorizability vanishes and 
suddenly reappears in the G-AF phase, not as a transi- 
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tion effect but rather as a robust feature vanishing only 
for high values of E z . We argue that this is related with 
surprisingly rigid interplane spin-spin correlations which 
should, if we think in spin-orbital factorizable way, decay 
quickly as t c approaches 1/2. Following this "factoriz- 
able reasoning" we could also expect stability of the C- 
AF phase for higher values of rj, above the G-AF phase. 
These effects are absent in our results, showing that intu- 
ition suggesting spin-orbital factorization can be mislead- 
ing even when considering such a simple isotropic orbital 
configuration. 

For higher values of Hund's exchange r) frustration in- 
creases when AF exchange interactions compete with FM 
ones, and as a result the most exotic phases with explicit 
on-site spin-orbital entanglement arise; two of them, the 
ESO and EPVB phase are neighboring and placed in be- 
tween the VBz and PVB ones, and they become degen- 
erate with both of them at the multicritical point where 
four phases meet (see Fig. [6]). This situation follows 
from the fact that singlet phases are more susceptible to 
ferromagnetism favored by high r\ than the G-AF phase 
is, which turned out to be surprisingly robust. The ESO 
phase is also a singlet phase similar to the VBz one but 
with spin singlets and orbital order gradually suppressed 
under increasing Hund's exchange r\. At the same time 
spin-orbital order stays almost constant and spin-orbital 
entanglement grows. 

Further increase of r\ always leads to the A-AF phase 
throughout a discontinuous transition accompanied by 
an abrupt drop of spin-orbital entanglement. Above rj ~ 
0.2 the ESO phase is completely immersed in the A-AF 
one and ends up with a single bicritical point. If we 
come back below r\ ss 0.2 then the ESO changes smoothly 
into the EPVB phase, being an entangled precursor of 
the PVB order, meaning that the non-uniform orbital 
order and in-plane singlets are formed and spin-orbital 
entanglement drops. On the other hand, this phase can 
be also seen as an extension of the A-AF into fully AF 
sector because the EPVB phase has long-range magnetic 
order, being however strongly non-uniform (see Fig. 14 1. 

Finally, we would like to remark that experimental 
phase diagrams of strongly correlated transition metal 
oxides are one of the challenging directions of recent re- 
search. Systematic trends observed for the onset of the 
magnetic and orbital order in the -RVO3 perovskites have 
been successfully explained by the competing interactions 
in presence of spin-orbital entanglement In contrast, 
the theory could not explain exceptionally detailed in- 
formation on the phase diagram of the i?Mn03 mangan- 
ites which accumulated due to impressive experimental 
workP^l The present K3CU2F7 bilayer system is some- 
what similar to La2_2a:Sri + 2 2 ;Mn207 bilayer manganites 
with very rich phase diagrams and competition between 
phases with different types of long-range order in doped 
systemsP^Such phases are generic in transition metal ox- 
ides and were also reproduced in models of bilayer man- 
ganites which have to include in a dditio n superexchange 
interaction between core ti g spin d 35 * 3 ^ that suppresses 



spin-orbital fluctuations and entanglement in the e g sub- 
system. In contrast, K3CU2F7 bilayer is rather unique as 
the only electronic interactions arise here due to entan- 
gled spin-orbital superexchange. They explain the origin 
of the VBz phase observed 3 ^ in K3CU2F7 but not found 
in bilayer manganites, and provide a unique opportunity 
of investigating whether signatures of spin-orbital entan- 
glement could be identified in future experiments. 

Summarizing, the presented analysis demonstrates 
that spin-orbital entanglement plays a crucial role in 
complete understanding of the phase diagram of the bi- 
layer spin-orbital d 9 model. By introducing additional 
spin-orbital order parameter independent of spin and or- 
bital mean fields we obtained phases with spin disorder 
in highly frustrated regime of parameters. The example 
of the entangled ESO and EPVB phases shows that joint 
spin-orbital order can be at least as strong as the other 
two (spin or orbital) types of order, or may even persist 
as the only symmetry breaking field when the remain- 
ing ones vanish. We argue that the cluster method wc 
used here is sufficiently realistic to investigate the phase 
diagram of the 3D spin-orbital d 9 model, and could be 
applied to other spin-orbital superexchange models ade- 
quate for undoped transition metal oxides. 
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Appendix: Solution of the mean-field equations 

Here we present briefly the solution of the self- 
consistency Eqs. ( 2.22 ) and ( 2.23 ) obtained in the single- 



site MF approximation. It is obtained as follows: assum- 
ing a—b symmetry of the system, i.e., putting x a = X b 
and £ a = £ b , we derive t a and t c from Eq. (2.21) as 
functions of a and /3, 
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Now we introduce a parametrization 

a — Asindi, ,3 = Acos< 
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t c one finds immediately sin</> depending on A, 



and use the self-consistency Eqs. (2.22) and (2.23). From 

(A.7) 



85A + I 



Comparing Eqs. (2.22) and (A.l) for t a one gets: 
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After inserting sin0 into Eq. (A.8) we obtain a surpris- 
ingly simple result for cos <j>: 



892 



0. 



(A.9) 



This leads to two classes of solutions of self-consistency 
Eqs. H2.22H and p.23); (i) either cos0 = 0, or (ii) A = 



3/8(?2 and cos <fi ^ 0. The first option implies sin</> = ±1 
and leads to two uniform orbital configurations with t c 



and t a = t" = -t c /2. Furthermore, using Eq. (A.7) 
we can calculate A and find the borders of these uniform 
phases demanding A > 0. 

The second option, i.e., cos(f> 7^ 0, implies AO-type of 
order with: 



tc 
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and with phase borders defined by the condition: 2\t c \ < 
1. The phase borders given here set the maximal range of 
the phase under consideration and cannot be treated as 
the lines of phase transitions shown in the phase diagram; 
the latter lines are determined by comparing the ground 
state energies E calculated form Eq. (2.24). 
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